Geospatial Carpentry | Session 3: Working with raster data
Claudiu Forgaci
2022-11-21
Setup
# Install packages
# library(remotes)
# remotes::install_github('ropensci/osmdata')
# Data manipulation
# library(tibble)
library(dplyr)
# library(reshape)
# Data visualisation
library(ggplot2)
# Working with spatial data
# library(sf)
# library(leaflet)
library(raster)
library(rgdal)
# library(proj4)
# library(scales)
# Working with OSM data
# library(osmdata)
# library(osmextract)Intro to raster data (30 + 20 minutes)
View Raster File Attirbutes
GDALinfo("data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif")## rows 1367
## columns 1697
## bands 1
## lower left origin.x 731453
## lower left origin.y 4712471
## res.x 1
## res.y 1
## ysign -1
## oblique.x 0
## oblique.y 0
## driver GTiff
## projection +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## file data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif
## apparent band summary:
## GDType hasNoDataValue NoDataValue blockSize1 blockSize2
## 1 Float64 TRUE -9999 1 1697
## apparent band statistics:
## Bmin Bmax Bmean Bsd
## 1 305.07 416.07 359.8531 17.83169
## Metadata:
## AREA_OR_POINT=AreaHARV_dsmCrop_info <- capture.output(
GDALinfo("data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif")
)Open a Raster in R
We are going to work with a Digital Surface Model (DSM) which is in the GeoTIFF format.
DSM_HARV <-
raster("data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif")
DSM_HARV## class : RasterLayer
## dimensions : 1367, 1697, 2319799 (nrow, ncol, ncell)
## resolution : 1, 1 (x, y)
## extent : 731453, 733150, 4712471, 4713838 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : HARV_dsmCrop.tif
## names : HARV_dsmCrop
## values : 305.07, 416.07 (min, max)summary(DSM_HARV)## HARV_dsmCrop
## Min. 305.5500
## 1st Qu. 345.6500
## Median 359.6450
## 3rd Qu. 374.2825
## Max. 413.9000
## NA's 0.0000summary(DSM_HARV, maxsamp = ncell(DSM_HARV))## HARV_dsmCrop
## Min. 305.07
## 1st Qu. 345.59
## Median 359.67
## 3rd Qu. 374.28
## Max. 416.07
## NA's 0.00DSM_HARV_df <- as.data.frame(DSM_HARV, xy = TRUE)str(DSM_HARV_df)## 'data.frame': 2319799 obs. of 3 variables:
## $ x : num 731454 731454 731456 731456 731458 ...
## $ y : num 4713838 4713838 4713838 4713838 4713838 ...
## $ HARV_dsmCrop: num 409 408 407 407 409 ...ggplot() +
geom_raster(data = DSM_HARV_df , aes(x = x, y = y, fill = HARV_dsmCrop)) +
scale_fill_viridis_c() +
coord_quickmap()
plot(DSM_HARV)
View Raster Coordinate Reference System (CRS) in R
But what units are these?
# CRS have been explained on Day 1
crs(DSM_HARV)## Coordinate Reference System:
## Deprecated Proj.4 representation:
## +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## WKT2 2019 representation:
## PROJCRS["WGS 84 / UTM zone 18N",
## BASEGEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]],
## CONVERSION["UTM zone 18N",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",-75,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["Engineering survey, topographic mapping."],
## AREA["Between 78°W and 72°W, northern hemisphere between equator and 84°N, onshore and offshore. Bahamas. Canada - Nunavut; Ontario; Quebec. Colombia. Cuba. Ecuador. Greenland. Haiti. Jamica. Panama. Turks and Caicos Islands. United States (USA). Venezuela."],
## BBOX[0,-78,84,-72]],
## ID["EPSG",32618]]Calculate the Min and Max value
minValue(DSM_HARV)## [1] 305.07maxValue(DSM_HARV)## [1] 416.07If Min and Max values haven’t been calculated, you can set them with
the raster::setMinMax() function.
# Maybe skip this
DSM_HARV <- raster::setMinMax(DSM_HARV)Raster bands
To see how many bands a raster dataset has, use the
raster::nlayers() function.
nlayers(DSM_HARV)## [1] 1This dataset has only 1 band. We will discuss multi-band raster data in a later episode.
Dealing with missing data
In raster data, pixels with no value are represented with
NoDataValue. That is usually the case with edge values
where data is missing from the rectangular region of the raster. Use
raster::GDALinfo() to check the metadata for missing data
values.
GDALinfo("data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif")## rows 1367
## columns 1697
## bands 1
## lower left origin.x 731453
## lower left origin.y 4712471
## res.x 1
## res.y 1
## ysign -1
## oblique.x 0
## oblique.y 0
## driver GTiff
## projection +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## file data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif
## apparent band summary:
## GDType hasNoDataValue NoDataValue blockSize1 blockSize2
## 1 Float64 TRUE -9999 1 1697
## apparent band statistics:
## Bmin Bmax Bmean Bsd
## 1 305.07 416.07 359.8531 17.83169
## Metadata:
## AREA_OR_POINT=AreaBad data values
Bad values are usually miscalculations that result in values out of a predefined range.
# This code chunk is hidden in the lesson materialCreating a histogram of raster values
A histogram can be used to inspect the distribution of raster values
visually. It can show if there are values above the max or below the min
of the expected range. We can create a histogram with the ggplot2
function geom_histogram().
ggplot() +
geom_histogram(data = DSM_HARV_df, aes(HARV_dsmCrop))
Adjust the level of desired detail by setting the number of bins.
ggplot() +
geom_histogram(data = DSM_HARV_df, aes(HARV_dsmCrop), bins = 40)
Challenge
GDALinfo("data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_DSMhill.tif")## rows 1367
## columns 1697
## bands 1
## lower left origin.x 731453
## lower left origin.y 4712471
## res.x 1
## res.y 1
## ysign -1
## oblique.x 0
## oblique.y 0
## driver GTiff
## projection +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## file data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_DSMhill.tif
## apparent band summary:
## GDType hasNoDataValue NoDataValue blockSize1 blockSize2
## 1 Float64 TRUE -9999 1 1697
## apparent band statistics:
## Bmin Bmax Bmean Bsd
## 1 -0.7136298 0.9999997 0.3125525 0.4812939
## Metadata:
## AREA_OR_POINT=AreaPlotting raster data (40 + 30 minutes - 10 minutes about plot formatting - 15 minutes details = 45 minutes)
DSM_HARV_df <- DSM_HARV_df %>%
mutate(fct_elevation = cut(HARV_dsmCrop, breaks = 3))
ggplot() +
geom_bar(data = DSM_HARV_df, aes(fct_elevation))
To see the cutoff values:
unique(DSM_HARV_df$fct_elevation)## [1] (379,416] (342,379] (305,342]
## Levels: (305,342] (342,379] (379,416]To show count number of pixels in each group:
DSM_HARV_df %>%
group_by(fct_elevation) %>%
count()## # A tibble: 3 × 2
## # Groups: fct_elevation [3]
## fct_elevation n
## <fct> <int>
## 1 (305,342] 418891
## 2 (342,379] 1530073
## 3 (379,416] 370835To customize cutoff values:
custom_bins <- c(300, 350, 400, 450)
DSM_HARV_df <- DSM_HARV_df %>%
mutate(fct_elevation_2 = cut(HARV_dsmCrop, breaks = custom_bins))
unique(DSM_HARV_df$fct_elevation_2)## [1] (400,450] (350,400] (300,350]
## Levels: (300,350] (350,400] (400,450]ggplot() +
geom_bar(data = DSM_HARV_df, aes(fct_elevation_2))
DSM_HARV_df %>%
group_by(fct_elevation_2) %>%
count()## # A tibble: 3 × 2
## # Groups: fct_elevation_2 [3]
## fct_elevation_2 n
## <fct> <int>
## 1 (300,350] 741815
## 2 (350,400] 1567316
## 3 (400,450] 10668ggplot() +
geom_raster(data = DSM_HARV_df , aes(x = x, y = y, fill = fct_elevation_2)) +
coord_quickmap()
Customising colors:
terrain.colors(3)## [1] "#00A600" "#ECB176" "#F2F2F2"ggplot() +
geom_raster(data = DSM_HARV_df , aes(x = x, y = y,
fill = fct_elevation_2)) +
scale_fill_manual(values = terrain.colors(3)) +
coord_quickmap()
my_col <- terrain.colors(3)
ggplot() +
geom_raster(data = DSM_HARV_df , aes(x = x, y = y,
fill = fct_elevation_2)) +
scale_fill_manual(values = my_col, name = "Elevation") +
coord_quickmap()
Challenge
DSM_HARV_df <- DSM_HARV_df %>%
mutate(fct_elevation_6 = cut(HARV_dsmCrop, breaks = 6))
unique(DSM_HARV_df$fct_elevation_6)## [1] (398,416] (379,398] (361,379] (342,361] (324,342] (305,324]
## Levels: (305,324] (324,342] (342,361] (361,379] (379,398] (398,416]my_col <- terrain.colors(6)
ggplot() +
geom_raster(data = DSM_HARV_df, aes(x = x, y = y,
fill = fct_elevation_6)) +
scale_fill_manual(values = my_col, name = "Elevation") +
coord_quickmap() +
xlab("X") +
ylab("Y") +
labs(title = "Elevation Classes of the Harvard Forest Digital Surface Model (DSM)")
Layering rasters
DSM_hill_HARV <- raster("data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_DSMhill.tif")
DSM_hill_HARV## class : RasterLayer
## dimensions : 1367, 1697, 2319799 (nrow, ncol, ncell)
## resolution : 1, 1 (x, y)
## extent : 731453, 733150, 4712471, 4713838 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : HARV_DSMhill.tif
## names : HARV_DSMhill
## values : -0.7136298, 0.9999997 (min, max)DSM_hill_HARV_df <- as.data.frame(DSM_hill_HARV, xy = TRUE)
str(DSM_hill_HARV_df)## 'data.frame': 2319799 obs. of 3 variables:
## $ x : num 731454 731454 731456 731456 731458 ...
## $ y : num 4713838 4713838 4713838 4713838 4713838 ...
## $ HARV_DSMhill: num NA NA NA NA NA NA NA NA NA NA ...ggplot() +
geom_raster(data = DSM_hill_HARV_df,
aes(x = x, y = y, alpha = HARV_DSMhill)) +
scale_alpha(range = c(0.15, 0.65), guide = "none") +
coord_quickmap()
ggplot() +
geom_raster(data = DSM_HARV_df ,
aes(x = x, y = y,
fill = HARV_dsmCrop)) +
geom_raster(data = DSM_hill_HARV_df,
aes(x = x, y = y,
alpha = HARV_DSMhill)) +
scale_fill_viridis_c() +
scale_alpha(range = c(0.15, 0.65), guide = "none") +
ggtitle("Elevation with hillshade") +
coord_quickmap()
Challenge
# CREATE DSM MAPS
# import DSM data
DSM_SJER <- raster("data/NEON-DS-Airborne-Remote-Sensing/SJER/DSM/SJER_dsmCrop.tif")
# convert to a df for plotting
DSM_SJER_df <- as.data.frame(DSM_SJER, xy = TRUE)
# import DSM hillshade
DSM_hill_SJER <- raster("data/NEON-DS-Airborne-Remote-Sensing/SJER/DSM/SJER_dsmHill.tif")
# convert to a df for plotting
DSM_hill_SJER_df <- as.data.frame(DSM_hill_SJER, xy = TRUE)
# Build Plot
ggplot() +
geom_raster(data = DSM_SJER_df ,
aes(x = x, y = y,
fill = SJER_dsmCrop,
alpha = 0.8)
) +
geom_raster(data = DSM_hill_SJER_df,
aes(x = x, y = y,
alpha = SJER_dsmHill)
) +
scale_fill_viridis_c() +
guides(fill = guide_colorbar()) +
scale_alpha(range = c(0.4, 0.7), guide = "none") +
# remove grey background and grid lines
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank()) +
xlab("UTM Easting Coordinate (m)") +
ylab("UTM Northing Coordinate (m)") +
ggtitle("DSM with Hillshade") +
coord_quickmap()
# CREATE DTM MAP
# import DTM
DTM_SJER <- raster("data/NEON-DS-Airborne-Remote-Sensing/SJER/DTM/SJER_dtmCrop.tif")
DTM_SJER_df <- as.data.frame(DTM_SJER, xy = TRUE)
# DTM Hillshade
DTM_hill_SJER <- raster("data/NEON-DS-Airborne-Remote-Sensing/SJER/DTM/SJER_dtmHill.tif")
DTM_hill_SJER_df <- as.data.frame(DTM_hill_SJER, xy = TRUE)
ggplot() +
geom_raster(data = DTM_SJER_df ,
aes(x = x, y = y,
fill = SJER_dtmCrop,
alpha = 2.0)
) +
geom_raster(data = DTM_hill_SJER_df,
aes(x = x, y = y,
alpha = SJER_dtmHill)
) +
scale_fill_viridis_c() +
guides(fill = guide_colorbar()) +
scale_alpha(range = c(0.4, 0.7), guide = "none") +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank()) +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank()) +
ggtitle("DTM with Hillshade") +
coord_quickmap()
Reproject raster data (40 + 20 minutes - 10 minutes instruction details - 20 minutes challenges = 30 minutes)
What happens when maps don’t line up? That is usually a sign that layers are in different coordinate reference systems (CRS).
DTM_HARV <- raster("data/NEON-DS-Airborne-Remote-Sensing/HARV/DTM/HARV_dtmCrop.tif")
DTM_hill_HARV <- raster("data/NEON-DS-Airborne-Remote-Sensing/HARV/DTM/HARV_DTMhill_WGS84.tif")DTM_HARV_df <- as.data.frame(DTM_HARV, xy = TRUE)
DTM_hill_HARV_df <- as.data.frame(DTM_hill_HARV, xy = TRUE)ggplot() +
geom_raster(data = DTM_HARV_df ,
aes(x = x, y = y,
fill = HARV_dtmCrop)) +
geom_raster(data = DTM_hill_HARV_df,
aes(x = x, y = y,
alpha = HARV_DTMhill_WGS84)) +
scale_fill_gradientn(name = "Elevation", colors = terrain.colors(10)) +
coord_quickmap()
ggplot() +
geom_raster(data = DTM_HARV_df,
aes(x = x, y = y,
fill = HARV_dtmCrop)) +
scale_fill_gradientn(name = "Elevation", colors = terrain.colors(10)) +
coord_quickmap()
ggplot() +
geom_raster(data = DTM_hill_HARV_df,
aes(x = x, y = y,
alpha = HARV_DTMhill_WGS84)) +
coord_quickmap()
### Challenge
crs(DTM_HARV)## Coordinate Reference System:
## Deprecated Proj.4 representation:
## +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## WKT2 2019 representation:
## PROJCRS["WGS 84 / UTM zone 18N",
## BASEGEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]],
## CONVERSION["UTM zone 18N",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",-75,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["Engineering survey, topographic mapping."],
## AREA["Between 78°W and 72°W, northern hemisphere between equator and 84°N, onshore and offshore. Bahamas. Canada - Nunavut; Ontario; Quebec. Colombia. Cuba. Ecuador. Greenland. Haiti. Jamica. Panama. Turks and Caicos Islands. United States (USA). Venezuela."],
## BBOX[0,-78,84,-72]],
## ID["EPSG",32618]]crs(DTM_hill_HARV)## Coordinate Reference System:
## Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs
## WKT2 2019 representation:
## GEOGCRS["WGS 84 (with axis order normalized for visualization)",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]],
## REMARK["Axis order reversed compared to EPSG:4326"]]Reproject rasters
DTM_hill_UTMZ18N_HARV <- projectRaster(DTM_hill_HARV,
crs = crs(DTM_HARV))crs(DTM_hill_UTMZ18N_HARV)## Coordinate Reference System:
## Deprecated Proj.4 representation:
## +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## WKT2 2019 representation:
## PROJCRS["WGS 84 / UTM zone 18N",
## BASEGEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]],
## CONVERSION["UTM zone 18N",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",-75,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["Engineering survey, topographic mapping."],
## AREA["Between 78°W and 72°W, northern hemisphere between equator and 84°N, onshore and offshore. Bahamas. Canada - Nunavut; Ontario; Quebec. Colombia. Cuba. Ecuador. Greenland. Haiti. Jamica. Panama. Turks and Caicos Islands. United States (USA). Venezuela."],
## BBOX[0,-78,84,-72]],
## ID["EPSG",32618]]crs(DTM_hill_HARV)## Coordinate Reference System:
## Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs
## WKT2 2019 representation:
## GEOGCRS["WGS 84 (with axis order normalized for visualization)",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]],
## REMARK["Axis order reversed compared to EPSG:4326"]]extent(DTM_hill_UTMZ18N_HARV)## class : Extent
## xmin : 731397.3
## xmax : 733205.3
## ymin : 4712403
## ymax : 4713907extent(DTM_hill_HARV)## class : Extent
## xmin : -72.18192
## xmax : -72.16061
## ymin : 42.52941
## ymax : 42.54234Challenge
Dealing with raster resolution
res(DTM_hill_UTMZ18N_HARV)## [1] 1.000 0.998res(DTM_HARV)## [1] 1 1DTM_hill_UTMZ18N_HARV <- projectRaster(DTM_hill_HARV,
crs = crs(DTM_HARV),
res = res(DTM_HARV))res(DTM_hill_UTMZ18N_HARV)## [1] 1 1res(DTM_HARV)## [1] 1 1DTM_hill_HARV_2_df <- as.data.frame(DTM_hill_UTMZ18N_HARV, xy = TRUE)ggplot() +
geom_raster(data = DTM_HARV_df ,
aes(x = x, y = y,
fill = HARV_dtmCrop)) +
geom_raster(data = DTM_hill_HARV_2_df,
aes(x = x, y = y,
alpha = HARV_DTMhill_WGS84)) +
scale_fill_gradientn(name = "Elevation", colors = terrain.colors(10)) +
coord_quickmap()
Challenge
# import DSM
DSM_SJER <- raster("data/NEON-DS-Airborne-Remote-Sensing/SJER/DSM/SJER_dsmCrop.tif")
# import DSM hillshade
DSM_hill_SJER_WGS <-
raster("data/NEON-DS-Airborne-Remote-Sensing/SJER/DSM/SJER_DSMhill_WGS84.tif")
# reproject raster
DTM_hill_UTMZ18N_SJER <- projectRaster(DSM_hill_SJER_WGS,
crs = crs(DSM_SJER),
res = 1)
# convert to data.frames
DSM_SJER_df <- as.data.frame(DSM_SJER, xy = TRUE)
DSM_hill_SJER_df <- as.data.frame(DTM_hill_UTMZ18N_SJER, xy = TRUE)
ggplot() +
geom_raster(data = DSM_hill_SJER_df,
aes(x = x, y = y,
alpha = SJER_DSMhill_WGS84)
) +
geom_raster(data = DSM_SJER_df,
aes(x = x, y = y,
fill = SJER_dsmCrop,
alpha=0.8)
) +
scale_fill_gradientn(name = "Elevation", colors = terrain.colors(10)) +
coord_quickmap()
Raster calculations (40 + 20 minutes - 15 minutes raster math - 15 minutes challenges = 30 minutes)
Raster calculations in R
GDALinfo("data/NEON-DS-Airborne-Remote-Sensing/HARV/DTM/HARV_dtmCrop.tif")## rows 1367
## columns 1697
## bands 1
## lower left origin.x 731453
## lower left origin.y 4712471
## res.x 1
## res.y 1
## ysign -1
## oblique.x 0
## oblique.y 0
## driver GTiff
## projection +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## file data/NEON-DS-Airborne-Remote-Sensing/HARV/DTM/HARV_dtmCrop.tif
## apparent band summary:
## GDType hasNoDataValue NoDataValue blockSize1 blockSize2
## 1 Float64 TRUE -9999 1 1697
## apparent band statistics:
## Bmin Bmax Bmean Bsd
## 1 304.56 389.82 344.8979 15.86147
## Metadata:
## AREA_OR_POINT=AreaGDALinfo("data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif")## rows 1367
## columns 1697
## bands 1
## lower left origin.x 731453
## lower left origin.y 4712471
## res.x 1
## res.y 1
## ysign -1
## oblique.x 0
## oblique.y 0
## driver GTiff
## projection +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## file data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif
## apparent band summary:
## GDType hasNoDataValue NoDataValue blockSize1 blockSize2
## 1 Float64 TRUE -9999 1 1697
## apparent band statistics:
## Bmin Bmax Bmean Bsd
## 1 305.07 416.07 359.8531 17.83169
## Metadata:
## AREA_OR_POINT=Area ggplot() +
geom_raster(data = DTM_HARV_df ,
aes(x = x, y = y, fill = HARV_dtmCrop)) +
scale_fill_gradientn(name = "Elevation", colors = terrain.colors(10)) +
coord_quickmap()
ggplot() +
geom_raster(data = DSM_HARV_df ,
aes(x = x, y = y, fill = HARV_dsmCrop)) +
scale_fill_gradientn(name = "Elevation", colors = terrain.colors(10)) +
coord_quickmap()
Raster math and Canopy Height Models (DROP THIS PART)
CHM_HARV <- DSM_HARV - DTM_HARV
CHM_HARV_df <- as.data.frame(CHM_HARV, xy = TRUE) ggplot() +
geom_raster(data = CHM_HARV_df ,
aes(x = x, y = y, fill = layer)) +
scale_fill_gradientn(name = "Canopy Height", colors = terrain.colors(10)) +
coord_quickmap()
ggplot(CHM_HARV_df) +
geom_histogram(aes(layer))
Challenge
custom_bins <- c(0, 10, 20, 30, 40)
CHM_HARV_df <- CHM_HARV_df %>%
mutate(canopy_discrete = cut(layer, breaks = custom_bins))
ggplot() +
geom_raster(data = CHM_HARV_df , aes(x = x, y = y,
fill = canopy_discrete)) +
scale_fill_manual(values = terrain.colors(4)) +
coord_quickmap()
Efficient raster calculations: overlay function
CHM_ov_HARV <- overlay(DSM_HARV,
DTM_HARV,
fun = function(r1, r2) { return( r1 - r2) })CHM_ov_HARV_df <- as.data.frame(CHM_ov_HARV, xy = TRUE) ggplot() +
geom_raster(data = CHM_ov_HARV_df,
aes(x = x, y = y, fill = layer)) +
scale_fill_gradientn(name = "Canopy Height", colors = terrain.colors(10)) +
coord_quickmap()
Export a GeoTIFF
writeRaster(CHM_ov_HARV, "CHM_HARV.tiff",
format="GTiff",
overwrite=TRUE,
NAflag=-9999)Challenge
CHM_ov_SJER <- overlay(DSM_SJER,
DTM_SJER,
fun = function(r1, r2){ return(r1 - r2) })CHM_ov_SJER_df <- as.data.frame(CHM_ov_SJER, xy = TRUE)ggplot(CHM_ov_SJER_df) +
geom_histogram(aes(layer))
ggplot() +
geom_raster(data = CHM_ov_SJER_df,
aes(x = x, y = y,
fill = layer)) +
scale_fill_gradientn(name = "Canopy Height",
colors = terrain.colors(10)) +
coord_quickmap()
writeRaster(CHM_ov_SJER, "chm_ov_SJER.tiff",
format = "GTiff",
overwrite = TRUE,
NAflag = -9999)ggplot(CHM_HARV_df) +
geom_histogram(aes(layer))
ggplot(CHM_ov_SJER_df) +
geom_histogram(aes(layer))
Work with multi-band rasters (40 + 20 minutes - 15 minutes details - 15 minutes challenges = 30 minutes)
##Getting Started with Multi-Band Data in R
The raster() function only reads in the first band, in
this case the red band of an RGB raster
RGB_band1_HARV <- raster("data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_RGB_Ortho.tif")RGB_band1_HARV_df <- as.data.frame(RGB_band1_HARV, xy = TRUE)ggplot() +
geom_raster(data = RGB_band1_HARV_df,
aes(x = x, y = y, alpha = HARV_RGB_Ortho)) +
coord_quickmap()
RGB_band1_HARV## class : RasterLayer
## band : 1 (of 3 bands)
## dimensions : 2317, 3073, 7120141 (nrow, ncol, ncell)
## resolution : 0.25, 0.25 (x, y)
## extent : 731998.5, 732766.8, 4712956, 4713536 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : HARV_RGB_Ortho.tif
## names : HARV_RGB_Ortho
## values : 0, 255 (min, max)nbands(RGB_band1_HARV)## [1] 3To import the green band:
RGB_band2_HARV <- raster("data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_RGB_Ortho.tif", band = 2)RGB_band2_HARV_df <- as.data.frame(RGB_band2_HARV, xy = TRUE)ggplot() +
geom_raster(data = RGB_band2_HARV_df,
aes(x = x, y = y, alpha = HARV_RGB_Ortho)) +
coord_equal()
Raster stacks
The stack() function brings in all bands
RGB_stack_HARV <- stack("data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_RGB_Ortho.tif")RGB_stack_HARV## class : RasterStack
## dimensions : 2317, 3073, 7120141, 3 (nrow, ncol, ncell, nlayers)
## resolution : 0.25, 0.25 (x, y)
## extent : 731998.5, 732766.8, 4712956, 4713536 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## names : HARV_RGB_Ortho.1, HARV_RGB_Ortho.2, HARV_RGB_Ortho.3
## min values : 0, 0, 0
## max values : 255, 255, 255RGB_stack_HARV@layers## [[1]]
## class : RasterLayer
## band : 1 (of 3 bands)
## dimensions : 2317, 3073, 7120141 (nrow, ncol, ncell)
## resolution : 0.25, 0.25 (x, y)
## extent : 731998.5, 732766.8, 4712956, 4713536 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : HARV_RGB_Ortho.tif
## names : HARV_RGB_Ortho.1
## values : 0, 255 (min, max)
##
##
## [[2]]
## class : RasterLayer
## band : 2 (of 3 bands)
## dimensions : 2317, 3073, 7120141 (nrow, ncol, ncell)
## resolution : 0.25, 0.25 (x, y)
## extent : 731998.5, 732766.8, 4712956, 4713536 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : HARV_RGB_Ortho.tif
## names : HARV_RGB_Ortho.2
## values : 0, 255 (min, max)
##
##
## [[3]]
## class : RasterLayer
## band : 3 (of 3 bands)
## dimensions : 2317, 3073, 7120141 (nrow, ncol, ncell)
## resolution : 0.25, 0.25 (x, y)
## extent : 731998.5, 732766.8, 4712956, 4713536 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : HARV_RGB_Ortho.tif
## names : HARV_RGB_Ortho.3
## values : 0, 255 (min, max)RGB_stack_HARV[[2]]## class : RasterLayer
## band : 2 (of 3 bands)
## dimensions : 2317, 3073, 7120141 (nrow, ncol, ncell)
## resolution : 0.25, 0.25 (x, y)
## extent : 731998.5, 732766.8, 4712956, 4713536 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : HARV_RGB_Ortho.tif
## names : HARV_RGB_Ortho.2
## values : 0, 255 (min, max)RGB_stack_HARV_df <- as.data.frame(RGB_stack_HARV, xy = TRUE)str(RGB_stack_HARV_df)## 'data.frame': 7120141 obs. of 5 variables:
## $ x : num 731999 731999 731999 731999 732000 ...
## $ y : num 4713535 4713535 4713535 4713535 4713535 ...
## $ HARV_RGB_Ortho.1: num 0 2 6 0 16 0 0 6 1 5 ...
## $ HARV_RGB_Ortho.2: num 1 0 9 0 5 0 4 2 1 0 ...
## $ HARV_RGB_Ortho.3: num 0 10 1 0 17 0 4 0 0 7 ...ggplot() +
geom_histogram(data = RGB_stack_HARV_df, aes(HARV_RGB_Ortho.1))
ggplot() +
geom_raster(data = RGB_stack_HARV_df,
aes(x = x, y = y, alpha = HARV_RGB_Ortho.2)) +
coord_quickmap()
Create a three-band image
plotRGB(RGB_stack_HARV,
r = 1, g = 2, b = 3)
Which is the same with terra
terra::plotRGB(RGB_stack_HARV,
r = 1, g = 2, b = 3)
plotRGB(RGB_stack_HARV,
r = 1, g = 2, b = 3,
scale = 800,
stretch = "lin")
plotRGB(RGB_stack_HARV,
r = 1, g = 2, b = 3,
scale = 800,
stretch = "hist")
Challenge
GDALinfo("data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_Ortho_wNA.tif")## rows 2317
## columns 3073
## bands 3
## lower left origin.x 731998.5
## lower left origin.y 4712956
## res.x 0.25
## res.y 0.25
## ysign -1
## oblique.x 0
## oblique.y 0
## driver GTiff
## projection +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## file data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_Ortho_wNA.tif
## apparent band summary:
## GDType hasNoDataValue NoDataValue blockSize1 blockSize2
## 1 Float64 TRUE -9999 1 3073
## 2 Float64 TRUE -9999 1 3073
## 3 Float64 TRUE -9999 1 3073
## apparent band statistics:
## Bmin Bmax Bmean Bsd
## 1 0 255 107.83651 30.01918
## 2 0 255 130.09595 32.00168
## 3 0 255 95.75979 16.57704
## Metadata:
## AREA_OR_POINT=AreaHARV_NA <- stack("data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_Ortho_wNA.tif")plotRGB(HARV_NA,
r = 1, g = 2, b = 3)
GDALinfo("data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_RGB_Ortho.tif")## rows 2317
## columns 3073
## bands 3
## lower left origin.x 731998.5
## lower left origin.y 4712956
## res.x 0.25
## res.y 0.25
## ysign -1
## oblique.x 0
## oblique.y 0
## driver GTiff
## projection +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## file data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_RGB_Ortho.tif
## apparent band summary:
## GDType hasNoDataValue NoDataValue blockSize1 blockSize2
## 1 Float64 TRUE -1.7e+308 1 3073
## 2 Float64 TRUE -1.7e+308 1 3073
## 3 Float64 TRUE -1.7e+308 1 3073
## apparent band statistics:
## Bmin Bmax Bmean Bsd
## 1 0 255 NaN NaN
## 2 0 255 NaN NaN
## 3 0 255 NaN NaN
## Metadata:
## AREA_OR_POINT=AreaRasterStack vs. RasterBrick
object.size(RGB_stack_HARV)## 51296 bytesRGB_brick_HARV <- brick(RGB_stack_HARV)
object.size(RGB_brick_HARV)## 170898912 bytesplotRGB(RGB_brick_HARV)
methods(class=class(RGB_stack_HARV))## [1] ! != [
## [4] [[ [[<- [<-
## [7] %in% == $
## [10] $<- addLayer adjacent
## [13] aggregate all.equal animate
## [16] approxNA area Arith
## [19] as.array as.character as.data.frame
## [22] as.integer as.list as.logical
## [25] as.matrix as.vector atan2
## [28] bbox blockSize boxplot
## [31] brick calc cellFromRowCol
## [34] cellFromRowColCombine cellFromXY cellStats
## [37] clamp click coerce
## [40] colFromCell colFromX colSums
## [43] Compare coordinates corLocal
## [46] couldBeLonLat cover crop
## [49] crosstab crs<- cut
## [52] cv density dim
## [55] dim<- disaggregate dropLayer
## [58] extend extent extract
## [61] flip freq getValues
## [64] getValuesBlock getValuesFocal hasValues
## [67] head hist image
## [70] init inMemory interpolate
## [73] intersect is.factor is.finite
## [76] is.infinite is.na is.nan
## [79] isLonLat KML labels
## [82] length levels levels<-
## [85] log Logic mask
## [88] match Math Math2
## [91] maxValue mean merge
## [94] metadata minValue modal
## [97] mosaic names names<-
## [100] ncell ncol ncol<-
## [103] nlayers nrow nrow<-
## [106] origin origin<- overlay
## [109] pairs persp plot
## [112] plotRGB predict print
## [115] proj4string proj4string<- quantile
## [118] raster rasterize ratify
## [121] readAll readStart readStop
## [124] reclassify rectify res
## [127] res<- resample rotate
## [130] rowColFromCell rowFromCell rowFromY
## [133] rowSums sampleRandom sampleRegular
## [136] scale select setMinMax
## [139] setValues shift show
## [142] spplot stack stackSelect
## [145] stretch subs subset
## [148] Summary summary t
## [151] tail text trim
## [154] unique unstack values
## [157] values<- weighted.mean which.max
## [160] which.min whiches.max whiches.min
## [163] wkt writeRaster xFromCell
## [166] xFromCol xmax xmax<-
## [169] xmin xmin<- xres
## [172] xyFromCell yFromCell yFromRow
## [175] ymax ymax<- ymin
## [178] ymin<- yres zonal
## [181] zoom
## see '?methods' for accessing help and source codemethods(class=class(RGB_stack_HARV[1]))## [1] [ [<- anyDuplicated as_tibble as.data.frame
## [6] as.raster bind boxplot brick cellFromXY
## [11] coerce coordinates determinant distance duplicated
## [16] edit extent extract head initialize
## [21] isSymmetric Math Math2 Ops raster
## [26] rasterize relist subset summary surfaceArea
## [31] tail trim unique values<- weighted.mean
## [36] writeValues
## see '?methods' for accessing help and source code